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Abstract 

We describe a method to simulate the dynamics of charged colloidal particles sus- 
pended in a liquid containing dissociated ions and salt ions. Regimes of prime current 
interest are those of large volume fraction of colloids, highly charged particles and 
low salt concentrations. A description which is tractable under these conditions is 
obtained by treating the small dissociated and salt ions as continuous fields, while 
keeping the colloidal macroions as discrete particles. For each spatial configuration 
of the macroions, the electrostatic potential arising from all charges in the system 
is determined by solving the nonlinear Poisson-Boltzmann equation. From the elec- 
trostatic potential, the forces acting on the macroions are calculated and used in a 
Brownian dynamics simulation to obtain the motion of the latter. The method is 
validated by comparison to known results in a parameter regime where the effective 
interaction between the macroions is of a pairwise Yukawa form. 
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1 Introduction 



Systems of charged particles of mesoscopic size suspended in a liquid are abun- 
dant in fields ranging from biology to chemical engineering [?]. Next to colloids, 
examples include proteins, polyelectrolytes, micelles etc. [?,?]. However, de- 
spite great interest from both fundamental and applied points of view, the 
phase behavior and transport properties of such systems are still only partly 
understood. 

The mesoscopic particles become charged because in a polar solvent like water 
small ions dissociate off of their surface. In addition to these dissociated ions 
the solution often also contains salt ions. In the following, no distinction will 
be made between dissociated ions and salt ions carrying the same charge. All 
the small ions are of molecular size and carry one or at most a few elementary 
charges. The mesoscopic particles, in contrast, have a much higher charge 
on them, hence, they are commonly referred to as macroions. Due to the 
electrostatic interaction, the oppositely charged counter-ions in the solution 
concentrate around the macroions while the like charged co-ions are depleted. 
Entropy opposes this tendency. As a result, an inhomogeneous distribution 
of small ions is established and part of the bare charge on the macroions 
is screened. The charges on the surface of the macroion together with the 
screening mobile charges in the solution are commonly referred to as the diffuse 
electric double layer [?,?]. 

When the double layers of two or more macroions overlap, the latter experience 
an effective force. Attempts to predict the macroscopic suspension properties 
mostly start from an assumed form of this interaction which is taken as pair- 
wise additive [?,?]. The validity of this assumption depends on the range of 
the double layer interaction which varies with the salt concentration and the 
charge on the macroions. For high salt concentration, the interaction range is 
small compared to the average distance between the macroions so that no more 
than two double layers overlap simultaneously. Thus, the double layer interac- 
tion is indeed pairwise additive and according to DLVO theory [?,?] described 
by a repulsive Yukawa potential. The effect of the macroionic charge in this 
case can be accounted for by renormalizing charge and size of the macroions 
[?,?,?]. At low salt concentrations, where the interaction range is several times 
the average distance between the macroions, it is very likely that the double 
layers of more than two macroions overlap at the same time. Hence, many- 
body effects become important and the description of suspension properties 
by a pair-potential becomes questionable. An explicit form for the additional 
many-body forces that arise, however, is not known at present. 

In the broad range of parameters where the form of the effective interaction 
is not known, a theoretical description must include the small ions explic- 
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itly. In the primitive model [?] both macroions and small ions are treated as 
spherical particles with different size and charge suspended in a continuous 
polarizable solvent. The feasibility of simulations using this model is limited by 
the number of particles involved: For highly charged macroions, many small 
ions are necessary to satisfy the constraint of overall charge neutrality and 
significant amounts of salt present in the solution further increase the number 
of particles. The most efficient methods currently allow simulation of only 80 
mesoscopic particles each carrying no more than 60 elementary charges and 
with no salt added to the solution [?]. In practice, the macroions carry up to 
several thousand elementary charges and salt concentrations between micro- 
moles and millimoles per liter are encountered while the volume fraction of 
the macroions is at most a few percent. Clearly, under these conditions prim- 
itive model calculations become intractable for reasonably large systems and 
a different approach is called for. 

A reduced description of the system is obtained by exploiting the fact that the 
small ions, having a much smaller mass, move on a much faster time scale than 
the macroions. This permits an adiabatic approximation, where partial equi- 
librium of the small ions for the instantaneous configuration of the macroions 
is assumed. Consequently, the small ions are described by a continuous density 
which minimizes an appropriate thermodynamic potential [?,?,?,?]. Using this 
density in the Poisson equation for the electrostatic potential results in a closed 
relation for the latter. When correlations between the small ions are neglected, 
this relation becomes the nonlinear Poisson-Boltzmann (PB) equation [?,?]. 
This mean-field description of the electrostatic effects has been shown to be a 
good approximation for the case of monovalent small ions and weak coupling 
between the small ions by comparison to a cell model Monte-Carlo simulation 
[?]. A more complete description of the electrostatics can be obtained by using 
more sophisticated free energy functionals [?,?]. 

Even on the PB level of description, the physics of colloidal suspensions 
presents a challenging problem so that further approximations are commonly 
invoked. At high volume fraction, the many-body problem posed by the sus- 
pension is reduced to an approximate one-body problem by using cell models 
[?,?]. However, the full nonhnear PB equation can be solved analytically only 
in very few cases [?]. Even for the extremely simple geometry of a spherical 
cell, no analytic solution is available. Therefore, its linearized version has of- 
ten been used as a substitute in deducing suspension properties [?,?]. By the 
combination of modern high-speed supercomputers and efficient numerical al- 
gorithms a treatment of the full nonlinear many-body problem is now within 
reach. 

In this paper, a method combining a PB field description of the small ions with 
a Brownian dynamics simulation of the macroions [?,?] is discussed. For each 
spatial configuration of the macroions, the electrostatic potential due to all 
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charges in the system is determined by solving the PB equation with bound- 
ary conditions determined by the positions of the macroions and the charges 
on them. The force acting on each of the macroions is then calculated by in- 
tegrating the stress tensor over a surface enclosing the respective macroion. 
These forces are finally used in a Brownian dynamics simulation to obtain 
the motion of the macroions. This method allows to calculate structural and 
thermodynamic properties of charge-stabilized colloidal suspensions with full 
account of the many-body interaction between the colloids mediated by the 
screening ionic fluid around them. 

To validate the method we here focus on a parameter range where the as- 
sumption of effectively pairwise additive interactions is expected to hold true. 
Specifically, we determine the effective pair-force for fixed arrangements of the 
macroions from the numerical solution of the PB equation and compare the 
results to linearized DLVO theory and predictions based on a cell model. This 
provides a stringent test of the crucial step in our simulation method, namely 
the calculation of the true forces acting on the macroions which serve as input 
for the Brownian dynamics part. The latter is then used to locate the melting 
point which for Yukawa systems is known from the classic work of Robbins, 
Kremer and Grest [?]. The excellent agreement for all of our comparisons in 
the regime of pairwise additive effective interactions confirms our method and 
implementation. Outside this parameter regime, many-body effects become 
important which are investigated elsewhere [?,?]. 

The paper is organized as follows: In Sect. [2] we define the system considered, 
collect the equations of motion, and discuss the different parameter regimes. 
The discretization of the PB equation and the method of solution are described 
in Sect. [3] while the force calculation and the Brownian dynamics algorithm 
are summarized in Sect. H] and Sect. [5l respectively. Sect. [6] presents results 
validating the method and illustrating its range of applicability. In Sect. [7] we 
finally discuss the results and point out directions for future development. 



2 System and Equations 

We consider N spherical macroions ( cf. Fig. [1]) in a cubic box of size L which 
is continued periodically in all three space directions to minimize finite size 
effects. The macroions have a radius a and and carry a charge —Ze which is 
distributed evenly over their surface. Two kinds of small ions with charge ±e 
may enter the system from a reservoir with ion pair concentration Cg. Thus, the 
number of small ions is not fixed but varies such that the chemical potential 
is constant everywhere in the suspension. In the absence of charges on the 
macroions {Z = 0), the average concentration of each species of small ions 
is Cg while otherwise their concentrations adjust so that charge neutrality is 
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Fig. 1. Illustration of a suspension with N = 3 macroions of charge —Ze suspended 
in a liquid of small ions with charges ±e. 

obeyed on average. Due to the Donnan effect, the total number concentration 
of small ions of any kind then is smaller than 2cs [?]. The presence of other, 
uncharged molecules in the solution is accounted for by a dielectric constant 
e. Finally, the whole system is coupled to a heat bath of temperature T. 

In addition to the macroion size a, there are several other important length 
scales which serve to distinguish different parameter regimes. Denoting the 
thermal energy as /3 = I/ZcbT, the Bjerrum length. 



gives the length below which the Coulomb interaction between two small ions 
dominates the thermal energy. Its value gives a measure of the importance of 
correlations between the small ions and, hence, the quality of the PB descrip- 
tion of the electrostatic problem. In case of monovalent small ions primitive 
model simulations and PB calculations for a spherical cell model have been 
compared in Ref. [?]. The deviations between both decrease with the value 
of the parameter A^/a and are already tiny at < 0.03a. For an aqueous 
solvent at room temperature, this corresponds to a > 24nm. Throughout our 
calculations we use an even smaller value Xs/a — 0.012 which corresponds to 
a particle size of a = dOnm as in recent experiments [?]. Hence, we are as- 
sured that PB theory furnishes an excellent description of the systems under 
investigation here. 

The width of the double layer and hence the range of the effective interaction 
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between the macroions is given by the Debye screening length k ^, where 



which emerges from a dimensional analysis of the PB equation. The volume 
available to each of the macroions, finally, yields a measure for the mean 
distance dm between two macroions as 



dm^N-y'L= [^y^, (3) 



where the volume fraction rj — has been introduced. For /t~^ <^ dm 

the interactions between the macroions are expected to be effectively pairwise 
additive. A sufficient condition that the potential ijj remains small everywhere 
is that in addition ZXb -C a. This may be used as a criterion ensuring that 
linearization of the PB equation will be valid. 



2.1 Electrostatic problem 



As discussed in the introduction, the Poisson-Boltzmann (PB) equation is 
obtained by combining the Poisson equation for the electrostatic potential with 
the charge densities of the two species of small ions resulting from minimization 
of a mean-field free energy functional J-'. The latter is composed of two terms 
representing the energy of all charges, fixed and mobile, in the electrostatic 
potential generated by them and the entropy for an ideal gas of small ions in 
contact with a particle reservoir 



/3^ = y infixed + n+ -n-) (4) 

G 

+ n,(.og(!^)-l) + „.(.og(^)-l)dV. 

Here, ±en±{r) are the number densities of the two species of small ions, 
■0(r) = f3e(f){r) is the scaled potential, and the integration volume G is the 
region outside the macroions. Minimization subject to the constraint of electro- 
neutrality gives 

n± = Csexp(:F^) . (5) 
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Together with the Poisson equation this leads to the PB equation which holds 
in the region G outside the macroions, i.e. 

V^'?/' = fi;^sinh(?/') for rEG, (6) 



where has been defined in equation [2J Inside the macroions, where there are 
no mobile charges, the electrostatic potential in principle is governed by Pois- 
son's equation and both regions are connected by a jump condition accounting 
for the charge distribution on the macroion surfaces. However, in aqueous so- 
lution, the dielectric constant inside the macroions is typically much smaller 
than that of the solvent outside and the matching condition between inside 
and outside simplifies to a von Neumann condition for the PB equation. De- 
noting the surface of the p-th macroion by dGp and its normal vector pointing 
into the solvent by Hp, this boundary condition is expressed as 

—ZX 

n„ ■ ^ijj = — for f e dGp . (7) 



Together with periodic boundary conditions on the sides of the box this fur- 
nishes a complete description of the electrostatic part of the problem from 
which other quantities can be derived. 



2.2 Force Calculation 



The forces on the particles needed for the Brownian dynamics simulation are 
calculated from the stress tensor, which is a sum of two parts, an osmotic and 
an electrostatic one. The osmotic part of the stress tensor 11 is proportional 
to the difference of the total number density of small ions from the reservoir 
concentration, i.e. 



pll = (n+ + n_- 2cs) 1 =2cs {cosh{ip) - 1) 1 
The electrostatic part of the stress tensor is 



PT = -— ((V^)'l -2VV^® VV^) . (9) 



Together we have 



Pt=P[u + t') (10) 
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((2/€2(cosh(^) - 1) + (VV^)') T-2VV'® V^) . (11) 



The force F"^ acting on the p-th macroion is obtained by integrating the normal 
component of the stress tensor Eq. flTT]) over the particle surface dG.p^ i.e. 

F^ = j T-fipdS . (12) 



It may be verified by straight forward calculation making use of the PB Eq. ([6]) 
that the divergence of the stress tensor vanishes. Therefore, the integration 
in Eq. (|T2|) needs not be taken over the particle surface; any other surface 
enclosing the macroion will give the same result. 



2.3 Colloid equation of motion 



Since the shape and charge distribution of the macroions are spherically sym- 
metric, only translational degrees of freedom have to be considered. Denoting 
the position of the p-th macroion hj Rp, p = 1 ... N , its equation of motion 
on the diffusive time scale follows from the force balance 

= F° + F/ + FS. (13) 



The electrostatic force F^ is calculated from the solution of the PB equation as 
described in the previous paragraph. Dissipative and stochastic forces model 
the heat bath. The dissipative forces F^ are given by 

= -C^p , (14) 



where ( = Q7ii]a is the Stokes friction coefficient of a macroion with radius a in 
the solvent of viscosity rj. The stochastic forces are related to the dissipative 
drag by the fluctuation dissipation theorem in order to ensure the correct 
equilibrium distribution. We have 



F; = ^2/cBTCep, (15) 



where T is the solvent temperature, is the Boltzmann constant and is 
an uncorrelated Gaussian white noise with zero mean and unit variance, i.e. 
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{Ut))=0 (16) 
2.4 Length and Time Scales 



To render the equations describing the colloid dynamics dimensionless it re- 
mains to introduce suitable scales to measure length and time. An obvious 
choice for the length scale is the macroion radius a. A physically sensible es- 
timate of the relevant time scale r is obtained as follows. An average force 
constant k for the interaction between the macroions may be defined in terms 
of the space-averaged Laplacian of the electrostatic energy as 



f3k=~Jv',pdV , (17) 

V 

where the volume of integration is the computational box. Using the PB equa- 
tion and the charge density of the small ions this becomes 



1 ^2 1 

pk= l{n+-n^)dV . 



3 2csV 

V 

Due to electro-neutrality, the integral must equal ZN, so that one finally finds 



Pk^^. (19) 

From this average force constant the time scale of the motion is obtained as 



r = JC/k. 



3 Spatial discretization and PB solution algorithm 



For large colloidal charge Z, the solution of the PB equation may develop 
rather steep boundary layers close to the particle surfaces while far away it 
varies only gradually (c/. Fig. [2]). To provide an accurate resolution of the 
boundary layers with a reasonable number of grid points, we follow Ref. [?] 
and define a spherical grid in a shell centered about each of the macroions 
overset on a Cartesian grid covering the simulation box. The solution of the PB 
equation on the domain G outside the macroions is then obtained in four steps 
(c/. Fig. [3]): First the PB equation is solved in each of the N spherical shells 
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Fig. 2. Potential landscape calculated for a random configuration of A*" = 5 
macroions in a cubic box, the size of which is given here in units of the grid spacing. 
The macroions are confined to the plane z = L/2 and the potential values in this 
plane are shown. The charge on each colloid is Z = 120 and the screening parameter 
is na = 0.5. 



fii, . . . , assuming given potential values at their outer edges. Interpolation 
of these spherical shell solutions to the Cartesian grid next yields boundary 
values for the interstitial region VLq which is the part of the domain G not 
contained in any of the spherical shells. With these boundary values the PB 
equation is then solved in the interstitial region. Finally, new boundary values 
on the outer edges of the spherical shells are found by interpolating back from 
the interstitial solution. These steps are repeated until a converged solution 
on the whole domain G is obtained. In the following we describe in turn the 
discretization of the PB equation on the Cartesian and spherical grids, the 
interpolation procedure between the grids, and the method used to solve the 
discrete equations. 
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L 



Fig. 3. Sketch of the computational box, a cube of side length L, with = 3 
macroions. The domain G outside the latter is partitioned into spherical shells Op 
around each macroion (crosshatched) and an interstitial region VLq comprising the 
rest of the simulation box (hatched). 

3.1 Discretized equations on the Cartesian background grid 



The cubic simulation box is covered by a Cartesian background grid with equal 
grid spacing A in the three space directions. The size of A is chosen so that 
it matches the radial step size at the outer edge of the spherical shells (c/. 
Sect. 13. 2p . Typical values resulting are A = 0.2. ..0.3a. 

Denoting by tpijk the value of the potential at the grid point with indices 
i,j, k in the x-, y-, and z— directions respectively, the PB Eq. IQ in Cartesian 
coordinates is discretized in a straight forward manner by the central difference 
approximations 



dipijk ^ 'ipi+ljk — i'i-ljk 

~dx 2A 

d'^i'ijk ^ ipi+ijk - ^ijjjjk + ipi-ijk 
dx^ ^ A2 

d'^i'ijk ^ "'Pi+lj+lk — i^i+lj-lk — ^i-lj+lk + i^i-lj-lk 

dxdy 4A2 ^ ' 
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and corresponding formulae for derivatives in the other coordinate directions. 
The mixed derivatives do not appear in the PB equation but are needed for 
the interpolation to the spherical regions (cf. Sect. 13.31) . The discretized PB 
equation in Cartesian coordinates then becomes 

Ipi+ljk + Ipi-ljk + -ipij+lk + -ipij-lk + ■ipijk+l + '^ijk-l - 'Olpijk 



Values of the potential at the overlapping grid points, i.e. those points of 
the Cartesian grid that lie in one of the spherical shells VLi are obtained by 
interpolation from the respective spherical grid as described in Sect. 13. 3[ These 
values are used as fixed boundary conditions for solving the PB Eq. f l21l) 
in the interstitial region VLq. The second set of boundary conditions for the 
Cartesian problem are the periodic boundary conditions on the surface of the 
computational box which are implemented in an obvious way. 

3.2 Discretized equations on the spherical grids around the particles 

In a local coordinate system with the origin at the center of a macroion the 
spherical shell around it is defined as the region between the macroion surface 
r = a and a spherical surface at some distance r = which is taken to be 
the same for all of the macroions. With the choice r^ < d/2, an overlap of 
a spherical region with another particle is avoided for practical purposes. To 
enhance the resolution near the particle surfaces, modified spherical coordi- 
nates are used, where the radial coordinate r is transformed to its inverse w 
according to [?] 



To fix the step sizes in the if-, 6'-, and ^-directions we first make a choice for 
so that the radial grid spacing at the particle surface, A^l r=a = 1/Ai„, is 
some small fraction of the particle radius, a typical value being Ar|r=a = 0.04a. 
Then Ag and A^ are adjusted so that at the equator of the coordinate system 
the grid points have approximately equal distances. Together with the value 
of Ts this also determines the number of grid points iV^,, Nq, and in the 
interior of the sliced spherical shell. The potential values at the grid points 
here are denoted as ifjijk = ipi'^'^w, j^e, kA^). Extra boundary points with 
indices i = 0, + 1 correspond to the outer edge of the spherical region and 
the particle surface, respectively. The points on the polar axis also require a 
special treatment and, hence, are labeled hj j = 0, Ne + 1. 

The PB Eq. ([6]) in the modified spherical coordinate system becomes 



(fi;A)^sinh (V'ijfc) • 



(21) 



r = 
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Using central differencing, its discretized form reads 



A ^ 2 



+ 



A? 



+ i^ijk+l . n A 

\Wi sm UjA^ 



1 



A, 
at 

A,, 



.WjAg 

sinh {^ijk) 



2wf tan 6'j Ae 
^ / A,„ 



+ -ipij-ik 
sin6'jA(^ 

\ 2n 



Wi sin 6*1 A 



A, 



A? 



2t(;?tan6'iAfi 



(24) 



On the polar axis of the spherical coordinate system, where ^ = or tt, 
the derivative with respect to is not defined and its coefficient in the PB 
Eq. (123|) diverges. A valid discretization at these points of the spherical grid is 
obtained from the integral version of the PB equation which by applying the 
Gauss statement becomes 



<fvil^-kdS = K^J sinh (tp) dV . (25) 

5 V 

A discrete relation for the grid points on the polar axis, tpijk with j = 0, Ng + 1, 
is obtained by considering a volume of integration bounded by the mid planes 
perpendicular to the grid lines connecting the point {i,j,k) to its nearest 
neighbors. For the points on the polar axis this is a frustum of pyramid the 
basis of which is formed by a regular polygon with A^,^ edges. Due to the 
smallness of the integration volume and its faces variations of ip in the volume 
and variations of on each face may be neglected. Approximating by the 
value at the grid point and by the value at the center of each face the 
following discrete expressions result: 



I kA,,, 



2 + 



2A„ 



wt 



sinh (ipio) 



2A^ / A^ 



TT \WiAe 



2 

^'ipiik 

k=l 



2A^ 



2A^ / A, 



2 



k=l 
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wf 



sinh {ipip^g+i) (26) 



from which the values of the potential at grid points on the polar axis are 
computed. Since at these points the polar angle is not defined, the third 
index k has been dropped. 

The boundary condition Eq. ([7]) on the particle surface translates into 

= ZXb , (27) 



dw 



w=l/a 

with the discretized form 



V'iYu.+ljfc — i^N^jk + ZXb^w , (28) 

where a one-sided difference approximation has been employed. The other set 
of boundary conditions are the values of the potential at the outermost surface 
r = Ts of the spherical region ipojk, which are determined by interpolation from 
the Cartesian grid. 



3.3 Interpolating between the coordinate systems 



After a pass on the Cartesian grid the values of the potential at the outer 
spherical boundary points with r = have to be related to the values just 
obtained for the Cartesian grid points. This is done in the following manner 
(c/. Fig. H]). For each point fgp on the edge of the spherical shell, first the 
closest point fca on the Cartesian grid is sought. Then all first and second order 
derivatives of the potential ip at that point are calculated according to Eq. (1201) 
using its 18 nearest and next-nearest neighbors. Introducing the components 
qa, a = x,y,z of the difference f^p - Vca, «-e. r^p - f^a = Y.a=x,y,z<laea the 
value at the spherical boundary point is then approximated as 



These values then serve as a constant-potential boundary condition for the 
solution of the problem in the spherical shell. 

Once the PB equation in the spherical shells has been solved, the values at the 
overlapping grid points, i.e. the Cartesian grid points falling inside a spherical 
region, have to be recalculated. To avoid the use of one-sided differences close 
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Fig. 4. At each outer boundary point fsp of the spherical shells, the value of the 

potential is defined by interpolation using the closest point fcaon the Cartesian grid 
(indicated by the arrow) and its nearest and ncxt-ncarcst neighbors (black squares). 
Note that in three dimensions this amounts to a total of 19 Cartesian grid points. 

to the edges of the spherical regions, this recalculation is done only for the 
Cartesian points lying inside a somewhat smaller sphere with radius = 
r^— Ar|r=rs/2, where A,.|,_rs/2 is the radial spacing at the edge of the spherical 
shell. For all of these points the potential values are interpolated in the same 
way as before in the interpolation from the Cartesian to the spherical grid. 
Note that even though points in the interior of the overlapping region never 
appear in the Cartesian grid calculation, their values may still be needed for 
the interpolations in case two or more spherical regions overlap. 

While in principle it may happen that a Cartesian grid point lies inside the 
particle itself, where the spherical grid is not defined, according to our expe- 
rience this hardly ever occurs in practice for suitably chosen parameters. To 
avoid uncontrolled abortion of the program, we nevertheless assign the average 
of all potential values on the particle surface to such points and monitor the 
frequency of their occurrence by a counter. 



3.4 Iterative Solution of the Discrete Equations 

On each of the grids the discretized PB equation is solved by a nonlinear 
Successive Over- Relaxation (SOR) method [?]. In each step of the iteration, 
first the residuum 
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Fig. 5. At the overlapping grid points (solid squares) of the Cartesian grid, the 
values of the potential are defined by interpolation from the spherical grid. 

TZijk = aipi+ijk + tnpi^ijk + cipij+ik + + eipijk+i + fipijk-i + g'^ijk 

-/isinh (V'ijfc) (30) 

is calculated at each interior grid point. Here a . . . h denote the coefficients 
in the discretized PB equation, Eq. fl2Tl) for the interstitial region Qq and 
Eq. ([23]) for the spherical shells fii, . . . , fijv- New values of the potential are 
used in the calculation as soon as they are available. Then the solution is 
updated according to 



ijjijkil + 1) = ^ijkil) - ' (31) 

where I is the iteration count. Finally, the boundary points are updated (c/. 
Eqs. (l26l) and (128|) ) except those where the value of the potential is fixed, of 
course) . 

For the over- relaxation parameter u we found that a value of = 1.5 was 
suitable for most cases, while sometimes it had to be reduced to = 1.3. 
The iteration is terminated when the maximum of the absolute value of the 
residuum for all grid points drops below some desired accuracy. More precisely 
the solution is considered converged when maXijk{\'JZijk/i^\) < 10~^. Note that 
it is not necessary to iterate to convergence in each of the repeatedly solved 
sub-domain problems for Qq, . . . ,Q]\r. Only in the last pass of this repetition, 
the full convergence must be ensured. In fact, we usually took only a single 
SOR step for each sub-domain problem. The necessary number of repetitions 
of the whole procedure then was a few hundred times for C(IOO^) Cartesian 
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grid points. 



4 Force Calculation 



As discussed in Sect. \2.2\ the electrostatic forces are obtained from an inte- 
gration of the normal component of the stress tensor over an arbitrary surface 
enclosing the particle. This integration is best performed in the spherical shell 
flp around each particle where the surface of integration is taken to be a sphere 
with radius r/ = l/wf. Introducing local base vectors 



/ sin 6 cos (p \ 




1 cos9cos(f)\ 




^ — sin \ 


sin 9 sin 


, ee = 


cos 9 sin (j) 


, = 


COS0 


V COS^ y 




\ — sin y 




^ ; 



(32) 



which are the same for the usual and our modified spherical coordinates, the 
normal vector and surface element are simply n = and dS = l/w"^. Thus, 
according to Eq. f|T2|) . the force becomes 



sin 9 



(33) 



The modified spherical components of the stress tensor are found from Eq. (|TT 
with the expression of the gradient 



2 dijj 4- ^ ^^p ^ ^ w dip , 



VV" = -w —e^ + w—ee + . 

ow o9 sm 9 



(34) 



The result is 



""^ sin 9 dw dcp 

A discrete approximation of Eqs. (15^ and fl5Bl) is obtained by using the mid- 
point rule for the integration and central differences for the derivatives of the 
potential. The best value for the radius of the integration surface = l/wj 
turned out to be somewhere in the middle of the spherical region. Close to the 
macroion surface, the variation of the potential is strongest so the solution of 
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the PB equation is less accurate there. Close to the outer edge of the spheri- 
cal region, on the other hand, some accuracy is lost due to the interpolation 
between the grids. 



5 Brownian dynamics 



The equation of motion Eq. ( |T3l) is advanced in time by a stochastic Euler 
algorithm [?,?] 



Rpit + h) = Rpit) + ^Fp{t) , (36) 

where the time step h in units of the time scale r defined in Sect. 12.41 is 
typically h = O.lr. The total force is given by 



4 = ^/ + l/^S, (37) 



where the last term is a time-discrete approximation to the white noise defined 
by Eq. (|T6l) . The Sp are random numbers with 

(Hp) = (38) 

drawn independently at each time step. Since we are interested only in calcu- 
lating averages from the trajectories, only the first two moments of the dis- 
tribution of Sp are important [?,?]. Hence, it is most convenient to generate 
these random numbers from a uniform distribution on the interval [—a/3, +v^] 
which has zero mean and unit variance. 



6 Validation of the method 



6. 1 Forces Between Two Isolated Macroions 



We first consider two isolated macroions, for which DLVO theory [?,?] predicts 
a repulsive double-layer interaction of Yukawa type with the potential 

f3U{r) = Uo— where Uo = i (39) 



r V 1 -|- Ka 
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Fig. 6. Dimensionless force / = PaF between two macroions placed in a large simu- 
lation box as function of their separation r/a (symbols) compared to the pair-force 
foLvo according to Eq. (HO]l predicted by linearized DLVO theory (solid lines). The 
simulation parameters are Z = 4 and Ka = 0.5 (circles), Z = 4 and Ka = 0.1 
(squares), and Z = 400 and na = 0.5 (triangles). Inset: the ratio f/fj^i^Yo ^ 
constant value of 1 (that is ^'o.{f / f^^^y^) = 0) when the simulation results agree 
with the DLVO force law, i.e. for small Z. Agreement up to a constant scale factor 
is still found for large Z as revealed by some other constant value ^'^{f / fuLvo) ^■ 

and the force 

PaF,,,,{r) = f/o(l + «:r)^e— . (40) 

These expressions are derived from the hnearized PB equation and hence 
expected to hold true if ZA^/a <^ 1. 

In Fig. [H] the force F calculated numerically from the solution to the full non- 
linear PB equation (symbols) is compared to the DLVO force law Eq. (1401) 
(solid lines). Good agreement between both is found for small macroionic 
charge Z = 4 (circles and squares). The slight deviations for the case with 
large screening length are due to the finite size of the computational box. For 
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large macroionic charge Z = 400 (triangles), however, there is a pronounced 
difference. By plotting the ratio F/F^^^^ it is seen that this difference lies 
entirely in the strength of the force and not in its dependence on the distance 
between the macroions. This suggests that even for large macroionic charge 
the force between two isolated macroions is of Yukawa form but with a renor- 
malized macroion charge Z^s instead of the bare charge Z {cf. Sect. 16. 3p . The 
validity of the linearization of the PB equation employed by DLVO theory 
can be checked directly by looking at the maximum value attained by the 
potential; the approximation sinh?/' ^ ip is valid for ip <^1. For the cases with 
Z = 4 the potential does not exceed 0.01 and 0.06, respectively for na = 0.5 
and Kfl = 0.1 so that the linearized PB equation is applicable. For the case 
with Z = 400, in contrast, the potential reaches a value of 5.5 and therefore a 
linearization of the PB equation cannot be justified. 

6.2 Effective Forces 

For the many-particle problem of a colloidal suspension it may be expected 
that the interactions effectively reduce to pairwise additive forces whenever 
the Debye length is small compared to the mean distance between the 
particles dm- The effective pair-forces are obtained for a fixed configuration 
of macroions in the following way [?]: Choosing two particles A and B, one 
first calculates the total force acting on particle B with particle A present. 
Then particle A is removed leaving all the other particles in place and one 
again calculates the total force F^ acting on particle B now with particle 
A removed. The force exerted by particle A on particle B then is just the 
difference between the two and varying the position of particle A results in 
the effective force curve: 

Fab{^ = Fl{r) - F° , f=RA-RB. (41) 

With this procedure, many-body interactions, if present, are folded into an 
effective pair interaction. If the true interactions in the system are pairwise 
additive, the resulting effective interaction is by construction identical to the 
true pairwise interaction potential. That is, it is independent of the direction 
of r and also independent of the arrangement of the surrounding particles (all 
particles but A and B). 

For macroions arranged in a FCC configuration, the resulting effective force 
calculated by numerical solution of the PB equation is shown in Fig. [71 The 
system parameters have been chosen to satisfy k,^^ <^ d where pairwise ad- 
ditivity is expected. The form of the interactions is again clearly seen to be 
of Yukawa type like in Eq. PUI) with some effective force parameters and 
Keff in place of the bare charge Z and the inverse Debye length n. 
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r/a 

Fig. 7. Dimensionless effective pair-force / = fiaFAB between two out oi N = 108 
macroions arranged in a FCC configuration. The volume fraction is r/ = 0.03 which 
gives a mean distance of = 5.2a. The bare charge of the macroions the screening 
parameter are Z = 3000 and Ka = 2.0. Note that for these parameters, Fab depends 
only on the particle separation r/a. Inset: Force multiplied by {r / a)^ / {1 + Hesr) and 
plotted logarithmically so that the Yukawa interaction appears as a straight line with 
slope —Kes- This scaling clearly shows that the interaction is Yukawa-like. From a 
fit, the effective force parameters are obtained as ZeS = 1080 and KeS = 1.98. 



To verify the pairwise additivity of the effective interaction this calculation 
has been repeated for the same parameters but taking different directions 
r = f/\f\ in a FCC configuration and also with the configuration of the sur- 
rounding particles changed to BCC. In all cases the same effective force law 
was found confirming that many-body interaction are indeed absent in this 
parameter regime. When the screening parameter Ka was reduced, however, a 
very different behavior results: The effective interaction becomes configuration 
dependent and exhibits a cut-off like feature at k,^^ ~ d. Such behavior is the 
manifestation of a many-body nature of the interactions as discussed in detail 
in Refs. [?,?]. 



21 



6.3 Charge Renormalization 
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Fig. 8. Effective force parameters Z^s and k^q for a Yukawa-type pair interaction 
as functions of the bare charge Z calculated for macroions in FCC and BCC config- 
urations (circles and squares respectively). The simulation parameters are = 108 
and 54 (for FCC and BCC configurations respectively), rj = 0.03, and na = 2.0. 
There is almost perfect agreement with the prediction of the spherical cell model 
[?] (dot-dashed lines). 

As shown in the previous section, for ^ the form of the effective 
interaction in a colloidal crystal is of the Yukawa type predicted by DLVO 
theory (c/. Eq. ( l40l) ). Except when ZXs/a is small, however, the bare charge Z 
and inverse screening length k have to be replaced by effective force parameters 
Zeff and Kefr to be determined by a fit to the calculated force curves. Repeating 
this procedure for different bare charges results in a charge renormalization 
curve Zcs as function of Z as shown in Fig. [8] (symbols) . For small bare charge, 
of course, Zes = Z as expected, but for large bare charge, Zes tends towards a 
constant value. The dependence of Kes on the bare charge is only very weak. 



The idea of charge renormalization has originally been proposed for a cell 
model [?], where only a single Wigner-Seitz cell of the colloidal crystal is 
considered. To facilitate solution of the PB equation, furthermore the com- 
plicated geometry of the true Wigner-Seitz cell is replaced by a sphere. The 
predictions of the spherical cell model for the effective parameters Z^s and Hes 
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Fig. 9. Electrostatic potential tp for N = 108 macroions arranged in a FCC con- 
figuration on a line corresponding to the 100 crystallographic direction with two 
particles located at x = 3.0 and x = 12.0 (symbols with solid lines to guide the 
eye). The system parameters are: volume fraction 77 = 0.023, screening parameter 
Ka = 0.3, and macroion charge Z = 21,41,61,81 from bottom to top curve. Inset: 
comparison to the potential calculated from the spherical cell model (dash-dotted 
lines). The agreement is good except near the cell edge where the spherical approx- 
imation becomes invalid. 



(dash-dotted lines in Fig. [8]) agree very well with the values obtained from the 
full numerical calculation (symbols). Again, there is no significant dependence 
on the crystal structure indicating the absence of many-body effects. 



A more direct comparison to the cell model can be made by looking at the 
electrostatic potential ij) itself. In Fig. [9] the numerical result for the potential 
along a line between two out of 108 macroions in an FCC configuration is 
shown for several values of the bare charge Z (symbols). The agreement to 
the spherical cell model result (dash-dotted line) is very good, except near the 
cell edge, where the spherical approximation leads to a higher value of the 
potential in the cell model. This is due to the fact that the half-width of the 
actual FCC Wigner-Seitz cell in the direction shown is 6 = 4.5 while the size 
of the spherical cell with equal volume is only R = 3.5. 
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Fig. 10. Mean square displacement {rms^) of N = 108 macroions as a function 
of time starting from an FCC configuration. Simulation parameters are: volume 
fraction rj = 0.03, screening parameter na = 1.0, and charges between Z = 100 and 
Z = 4000. The value of the scaled temperature lO^T as defined in Eq. ()44p is shown 
next to each curve. According to the Lindemann criterion, melting occurs whenever 
the root mean square displacement exceeds 19% of the mean distance between the 
macroions [?] (dashed line). 

6.4 Melting Point 

Finally, we take a look at a macroscopic property of the suspension obtained 
from a full dynamical calculation, namely the melting point. For systems of 
point-like Yukawa particles this quantity is known from the classic work of 
Robbins, Kremer, and Grest [?] . In this work, the melting line has been calcu- 
lated from a molecular dynamics simulation by using the Lindeman criterion 
stating that melting occurs whenever the root mean square displacement ex- 
ceeds the value of 19% of the mean distance between macroions. The result 
is 

f ^^(A) = 0.00246 + 0.000274A (42) 



in terms of the parameters 
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A = Kcsdm and (43) 

The former is simply the inverse effective screening length scaled with the 
mean particle distance. The latter represents a scaled temperature {kT in 
units of the Einstein phonon energy) which is related to the strength of the 
interaction Uo defined in Eq. (l39l) and, hence, the effective charge Zes- The 
function 6{X) depends on the lattice-type of the crystalline phase and is given 
in Tab. I of Ref . [?]. 

In Fig. [in] we show the mean square displacement, rms^, of the macroions as 
a function of time as calculated from the combined PB-Brownian dynamics 
method. By using the effective force parameters Zes and Keg from the previous 
section, the values of A and T can be calculated so that a comparison to 
Eq. becomes possible. The value of A = 8.1 is practically independent of 
the bare charge Z. The effective temperature, in contrast, depends strongly 
on Z with values shown next to each curve in Fig. [TOl The Lindemann value 
rms/dra = 19% is reached for T > 5.2 ■ 10~^, so that by the same criterion 
as in Ref. [?] the colloidal crystal melts around T = 5.1 ■ 10~^ ± 10%. The 
melting occurs quite close to the saturated value of the effective charge. The 
corresponding effective temperature T 4.5 ■ 10~^ cannot be exceeded for our 
system. We conclude that our estimate of the melting point is in reasonable 
agreement with the value 5.3 ■ 10"'^ obtained by Robbins et al. in Ref. [?] for 
Yukawa systems. 



7 Discussion 



The physics of charge stabilized colloidal suspensions poses challenging prob- 
lems resulting from the presence of many coupled degrees of freedom. A treat- 
ment has been possible so far only in limited regions of parameter space either 
due to the necessity of approximations in analytical treatments or because of 
practical limitations in numerical simulations. We have here described in detail 
a method which makes accessible the regime of large macroionic charge, low 
salt concentration, and relatively high colloid volume fraction, where many- 
body interactions between the macroions become important [?,?]. This has 
been achieved by combining a Poisson-Boltzmann field description of the small 
ions with a Brownian dynamics simulation of the charged colloidal particles. 
By describing the small ions in terms of a continuous density, the compu- 
tational effort becomes independent of their number. In this way, the main 
limitation, which precludes the application of primitive model simulations to 
the above parameter regime is circumvented. At the same time, in contrast to 
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analytical treatments like linearized DLVO theory or cell models, macroionic 
many-body effects are fully accounted for in our approach. 

In the present work, we have focused on the case where the effective inter- 
actions are still of a pairwise additive Yukawa form. By comparing our sim- 
ulation results to linearized DLVO theory and cell model calculations, which 
are applicable in this regime, the method has been thoroughly validated. Fur- 
thermore, its potential applications have been illustrated which reach beyond 
simple structural calculations [?,?] and also include the evaluation of thermo- 
dynamic properties, e.g. the phase behavior of charged colloidal suspensions. 
These applications are pursued in greater detail elsewhere [?]. 

Directions for further development from a physical viewpoint are twofold. 
First, on the Poisson-Boltzmann level of description correlations between the 
small ions are neglected, which limits the applicability to monovalent small 
ions and small coupling <^ a. The range of applicability could be extended 
by using more elaborate density functional theory similar to Refs. [?,?]. Sec- 
ondly, the calculation of rheological properties necessitates inclusion of hydro- 
dynamic effects. The framework of a field description of the solvent on which 
our approach is built, readily accommodates such an extension as well [?]. 

Desirable improvements of the computational efficiency include parallelization 
of the method and the implementation of more sophisticated solvers for the 
discrete equations. Several possibilities for the latter have been compared in 
Ref. [?] where a combination of inexact Newton and multigrid methods was 
found to be the most efficient. A parallelization strategy that naturally fits 
together with the use of the overset spherical grids is to use a domain decom- 
position approach [?] also for the Cartesian background grid. Developments 
along these lines are currently pursued. 
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